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Abstract 

A new order parameter approximation to Random Boolean Networks (RBN) 
is introduced, based on the concept of Boolean derivative. A statistical argu- 
ment involving an annealed approximation is used, allowing to measure the 
order parameter in terms of the statistical properties of a random matrix. 
Using the same formalism, a Lyapunov exponent is calculated, allowing to 
provide the onset of damage spreading through the network and how sen- 
sitive it is to minimal perturbations. Finally, the Lyapunov exponents are 
obtained by means of different approximations: through distance method and 
a discrete variant of the Wolf's method for continuous systems. 
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1 Introduction 



Random Boolean networks (RBN), also called Kauffman nets [1-2], were originally 
formulated as a model of genetic nets. The computational and mathematical prob- 
lems arising from continuous dynamical systems with a very high number of coupled 
non-linear equations lead to the introduction of RBN as an alternative approach. 
In this way, relevant statistical properties of RBN were derived [18]. This general 
analysis allowed to test several hypothesis concerning the large-scale organization of 
biological regulatory networks. 

Recent studies in the field try to obtain a natural bridge between discrete and 
more biologically sensible, continuous networks [3]. In this vein for instance, a 
characteristic quantitative measure associated to continuous dynamical systems is 
the Lyapunov exponent A [4]. This parameter measures the degree of instability 
of continuous dynamical systems, and allows us to characterize the transition to 
chaos by measuring the pace at which initial conditions tend to diverge as the 
system evolves. Although Lyapunov exponents have been also derived (or estimated) 
for discrete systems (as cellular automata [5-8]) there is, as far as we know, no 
study about this quantity in RBN and related systems. To have a way to estimate 
Lyapunov exponents for RBN's would thus establish a natural link between discrete 
and continuous systems. 

The paper is organized as follows. First the RBN formahsm and its order- 
disorder phase transition are introduced. Secondly, using the concept of the Boolean 
derivative [9] we propose a new order parameter for the RBN phase transition: the 
percent of I's in the Jacobian matrix that represents the Boolean derivative of the 
system. We then define a Lyapunov exponent for the RBN and compare our results 
with the distance method [11]. Finally, a second possible order parameter (the 
self-distance) is introduced. This will allow us, through a discrete analog of Wolf's 
method for continuous systems, to reobtain an expression for the Lyapunov exponent 
consistent with that previously found. 

2 Random Boolean Networks 

A RBN is a discrete system involving N units/ automat a with two possible states of 
a boolean variable {0, 1}. Each automaton is randomly connected with exactly K 
neighbors. The state of each unit is updated by means of a Boolean function, also 
randomly chosen from the set of all the Boolean functions with K binary arguments. 



2 




OR 


OR 


AND 


12 


3 




1 3 


2 




23 


1 


00 







00 







00 





01 


1 




1 


1 




1 





10 


1 




10 


1 




1 





1 1 


1 




1 1 


1 




1 1 


1 




100 



no — " on ^ni ) 

/ ^ 

101 



Figure 1: Example of RBN with N — 3 and K = 2. Left: we show the interac- 
tion graph and the corresponding rule-tables (Boolean functions) of each automata. 
Right: we show explicitly the transitions between global system states as flow dia- 
grams. 

Once the neighborhood and functions have been chosen they are fixed in time (i. e. a 
quenched set is used). The RBN exhibit a second order transition: for X < 2 a frozen 
(ordered) phase is observed, while for X > 2 a disordered phase sets in. A RBN is 
by definition a discrete (iV cells) deterministic system with a finite number of states 
({0, 1}), and therefore periodic patterns are expected after a maximum of 2^ steps. 
Thus, if we follow strictly the standard definition of low-dimensional deterministic 
chaos, chaotic behavior is not possible in these systems. Taking this into account, 
we will define chaotic behavior here through damage spreading [12-14]: a phase will 
be chaotic if damage spreading takes place, i. e. if changes caused by transient flips 
of a single unit propagate and grow until they reach a size comparable to that of 
the system. Thus our disordered phase will be called chaotic phase, analogously to 
continuous systems. 

Let us illustrate the RBN structure and dynamics with a simple example [2] 
wich will be used below. Given a system with = 3 automata with values: Xi,X2 
y X3 and connectivity K = 2, the net is wired by choosing the input neighbors as 
indicated in flgure 1. 
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The inputs have been represented in figure 1 as arrows that connect the automata 
and a simple directed graph is obtained. The state of the system at any time is an 
ordered array of bits, X = (xi, X2, 2:3), and the interactions between the automatas 
are described by Boolean functions. In this example the Boolean functions have 
been randomly sampled from the set of 2^^ = 16 possible functions with K = 2 
arguments: 

(1) For the automaton 1: /i(x2,X3) the function AND. 

(2) For the automaton 2: f2{xi,xz) the function OR. 

(3) For the automaton 3: f^{x\,X2) the function AND. 

These functions are represented in figure 1 (left) by means of rule-tables (all possi- 
ble inputs with their corresponding outputs) . The system is updated synchronously. 
A possible temporal succession of states will describe a trajectory (orbit) of the sys- 
tem. In figure 1 (right) we represent all possible trajectories of the system as a fiow 
diagram. 

In early studies on RBN phase transitions the critical point was estimated 
through numerical simulations [1,2]. The critical connectivity K = 2 gave the 
transition order-chaos. Later on this transition point was analytically obtained by 
means of the so-called Derrida's annealed approximation [6,11], also known as the 
distance method. Derrida developed a non-correlated (annealed) RBN model by 
randomizing the inputs and Boolean functions at each time step showing that, in 
the thermodynamic limit, the transition point is the same in the quenched and the 
annealed systems. This approach can be extended to a continuous (average) K- 
valued and biased RBN [15,16]. In biased RBN's the Boolean functions are chosen 
with a bias p, that is: the mean percentage of I's in the output is p. From two 
replicas with an initial normalized Hamming distance, d{t = 0), we can derive the 
equation for the evolution of d{t): 

d{t+l)^2p{l-p){l-[l-d{t)f} (1) 

At the frozen phase, the fixed point d* — is stable (i.e. the two initial con- 
figurations become identical as they evolve). In the chaotic phase however d* — 
is unstable, and two initially close configurations diverge to a finite distance. The 
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Figure 2: Phase space for RBN. The critical hne (continuous) is given by equation 
2. It separates the ordered and chaotic phases. For K = ?> (dashed hne), three 
particular examples are shown for N — and random initial conditions. Here 
Si — 1,0 are indicated as black and white squares, respectively. In all space-time 
diagrams, time runs from bottom to top. The three diagrams correspond to: p = 
0.60 (chaotic phase), p — 0.90 (ordered phase) and p — 0.79, at the transition line. 



critical curve on the parameter space (p, K) is: 

1 



2p{l-p) 



(2) 



In figure (2) the phase space is shown. On a constant connectivity line (here 
K = 3) three different runs of the system have been chosen for a RBN with K = 3 
and N = 50. For p = 1/2 equation (2) reduces to the standard RBN problem [2]. 
Fig. 3 shows with continuous lines the evolution of d{t) towards the theoretical fixed 
point d* (obtained from iteration of Eq. (1) with K — 3), p changing for each hne. 
The values chosen for p in this figure match the dots on the dashed K — 3 hne in Fig 
2. The squares represent the average result of 100 runs with two different replicas 
each. Here we used N — 10.000 automata, and the initial distance is d{0) — 0.5. 
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Figure 3: Continuous lines: dynamical evolution of the distance between two identi- 
cal replicas of RBN with initial distance of = .5 for K = ?> and different values of the 
bias p through the iteration of equation (1). Squares (points): numerical averages 
of the distance (auto-distance) for 100 experiments with RBN of size N — 10.000, 
K — 2> and the bias showed. 

Observe that d* acts as an order parameter. In fig. 4 the continuous line represents 
the stationary values obtained by iteration of equation (1) for changing p. Again, 
the squares represent the numerical values obtained by averaging over 100 runs with 
two different replicas of RBN each, with size N — 10.000, and d(0) = 0.50. 

In [17], Flyvbjerg defined a different order parameter: he defined the stable core 
at time t as the set of units that have reached stable values at time t (that is, remain 
unaltered in value for t' > t and are independent of the initial conditions). Let us 
define s{t) as the relative size of the stable core at time t, i.e. s{t)N is its absolute 
value. Then the asymptotic stable core size, s* = limt^oo s{t), is an order parameter 
for the order-chaos transition in RBN's. Flyvbjerg obtains an iterated equation for 
the stable core: 

s{t+l)=j:(^]s{tf-\l-s{t)yp, (3) 
1=0 \ ^ / 

where Pi is the probability that the Boolean function output be independent of a 
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certain number i of inputs. For the Boolean functions with bias p it yields pi — 
p'^' + (1 — pY' . By analyzing the stability of Eq. (3), he found out an identical 
transition curve than the one given by Eq. (2). In Fig. 4 (Short-dashed line) we 
represent this order parameter as 1 — s*. 



3 Boolean Derivatives in RBN 

We will now define a RBN in a more formal way. A RBN is a discrete dynamical 
system whose evolution is given by the iteration of a global mapping: 

F^,, : {0, 1}^ ^ {0, 1}^ (4) 

where Fk,p = (/i, /2, In), and with each /j being a Boolean function of K argu- 
ments and bias p (mean percentage of ones in the outputs) : 

fK,, : {0, 1}^ ^ {0, 1}, (5) 

A given configuration at time t, x(t), is updated sincrhonously, i.e.: 

x*+i = FkA^') (6) 

where each automata with £ {0, 1} is updated by mean of its corresponding 
Boolean function: 



X 



For a given Fx^pit) we define its x Jacobian matrix, F^p(t), as that whose 
elements are given by the Boolean derivatives at time t [9]: 

F' .(t) = ^f^^^^i^ = I Mxl,, ■■■,x/, ...,xU ® fi{xl^, ...,x/, ...,xU if Xj inputs Xi 
''^ dx/ lo otherwise 

(8) 

Here © is the exclusive OR (XOR) Boolean operation and Xj — Xj ® 1 (i.e., the 
binary complement of Xj ) . From the point of view of damage spreading [7] , we can 
see that j(t) = 1 if a flip in the input a;* at time t generates a change of x*"*"^ to 
x*"*"^ in step t + 1. In others words, the function spreads the damage. Otherwise, 
Flj{t) = 0, and no damage is spread. Note that F'i^p{t) depends on t because its 
value depends on the concrete configuration at time t. 

Continuing with our example (see Fig 1), let us suppose that, at t, the state of 
the system is x(i) = (1, 0, 1). If we compute the 3x3 Jacobian matrix F2 Q ^{t) its 
components will be: 

^1 = 
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^l'2 = /l(0,l)©/l(l,l) 

i^i',3 = /i(0,l)©/i(0,0) 

F2',i = /2(1,1) 0/2(0,1) 

F2% = /2(1,1)®/2(1,0) 
^3,1 = /3(l, 0)0/3(0,0) 
i^3,2 = /3(l,0)®/l(l,l) 
^3,3 = 

Thus the Jacobian matrix is: 

/O 1 0\ 
F'{t) = 

VI 0) 

Analogously to the continuous dynamical systems counterpart we take a config- 
uration state x(t) and a slight perturbation of it y(t). A perturbed configuration 
is a new configuration at a non-zero (but otherwise small) Hamming distance from 
the original one. In fact, it is possible to define the perturbation d{t) such that: 

y(t)=x(t)©d(t) (9) 

where the normahzed Hamming distance between x{t) and y{t) is 

1 ^ 

\mh^T.4 (10) 

i=l 

In our example, we have x(t) = (1, 0, 1) and we will take as the perturbed configu- 
ration y(t) = (0, 0, 1). Thus: 

y{t) = (0, 0, 1) = x(i) ® d{t) = (1, 0, 1) ® (1, 0, 0) = (1 ® 1, ® 0, 1 © 0) = (0, 0, 1) 
Note that we have the minimum possible perturbation 

\d{t) |=i(o + o + l) = l 
and that we can write the perturbation as 

d(i)=x(t)®y(i) (11) 



= 0® 1 = 1 
=0©0=0 
=10 1=0 

= 1©1 =0 
=100=1 
=101=0 
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Thus, in our example: 



(1, 0, 1) © (0, 0, 1) = (1 © 0, © 0, 1 © 1) = (1, 0, 0) 

Now we are ready to find the approximate evolution of the perturbation, as it is done 
in continuous systems. Using the equations (6), (9), (11), and using the Jacobian 
matrix to make a linear approximation, we have: 

d{t + 1) = y{t + 1) © x(i + 1) = F(y(t)) © F(x(i)) = 

= F(x(t) © d{t)) © F(x(t)) F'(x(t)) d{t) (12) 
where we define as: 

F'(x(t)) d{t) = e(F'(x(t)) . d{t)) (13) 

and the vector 0(x) = {Q{xi), Q{xi), ...,Q{xn)) having standard Heaviside func- 
tions as components: 

^ ^ ll otherwise ^ ^ 

Thus, in our example we have: 
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4 Jacobi matrix and a new order parameter 

Our aim is now to present a new order parameter based on the Jacobian matrix 
of Boolean derivatives, as previously defined. The dynamical equation (12) can 
be iterated in t to determine the evolution of the perturbation with two possible 
outcomes: as t — cxo either the initial perturbation at t = will tend to disappear, 
I d(cxo) 1^ 0, or it will reach a finite value. The behavior of the perturbation will 
be determined by the successive products of the Jacobian matrix. Thus, we define: 

M{t) = F'(x(0)) F'(x(l)) ... F'(x(i)) (15) 

If the number of I's in the Jacobian matrix is small the product in (15) will converge 
to a matrix M* formed only by zeros, and any initial perturbation will disappear. 
We will now attempt to construct an iterated equation for the evolution of the frac- 
tion of zeroes in the matrix M using a mean field approach. If our system has a 
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Figure 4: Continuous line: d* are assymptotic distances reached by iteration of 
equation (2) for different p with K = 3. Short-dashed hne: s* is the assymptotic 
unitary percent of elements of the stable core reached by iteration of equation (3) by 
changing p with K = 3. Long-dashed line: 1 — q* is the asymptotic unitary percent 
of ones in the matrix M* by iteration of equation (16). 



connectivity K and bias p, we substitute at each time step (mean field approxima- 
tion) the deterministic matrix F'(t) by a random matrix Q of the same form (at 
most K I's at each row). The probability for a randomly generated Boolean func- 
tion to have a Boolean derivative with value one F^^ = 1 is equal to the probability 
that a flip in its input Xi^ generates a change in its output Xi to Xi. We have two 
possibilities: the output has a value 1 and changes to 0, with probability p{l — p) 
and the symmetric case with probability {l—p)p. Thus the mean number of I's per 
row is 2p(l —p)K. At each time step, i -|- 1, we multiplie M(t) by a random matrix 
with a mean number 2p{l — p)K of I's per row, i.e. M(t + 1) = i7M(t). 

This approach is clearly analogous to the distance method [11], where the Boolean 
functions and inputs of the system vary at random at each time step. So, if we as- 
sume that at a time step t the matrix M(i) has a percentage qt of zeros, in the 
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Figure 5: Open circles: values of Lyapunov exponent define from equation (27) with 
K = ?> and different p variating. Filled circles: numerical estimation of the Lyapunov 
exponent by avering the expansion rate (18) calculate through the distance equation 
(1) with T = 10 in equation (19). 



thermodynamic limit, the g^+i will be: 



qt+i = lim 

Tv^oo 



^ 2p{l-p)K 
N 



(1 - <lt) 



N 



^-{l-qt)2p(l-p)K 



(16) 



Where 2p(l — p)K/N is the probability that flij = 1 and 1 — is the probability 
that Mjk = 1. By analyzing the stability of (16) around q* = 1 (i.e. all the elements 
of M* are zero), we find out the critical transition curve, given again by 



't+i 



dqt 



K2p{l-p) < 1 



;i7) 



q*=l 



in agreement with previous results [7,11,17]. In Fig. 4 we indicate (long-dashed 
hne) this order parameter as 1 — g* (the percentage of I's in the infinite product of 
Jacobians) for K — 3. 
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5 Lyapunov exponents and RBN 



In the previous section, we have introduced the Boolean derivative F'(x(t)), in 
analogy with the standard continuous counterpart. This operator was then used 
in order to define the discrete map (12) which gives us the time evolution of the 
perturbation d{t). Using this definition, an expansion rate of perturbations for RBN 
can be easily defined. The damage expansion rate will be [6]: 



d(t+ 1) 
I d{t) I 



(18) 



This allows us to define a Lyapunov exponent: 



A(T) = -Elogr7(i) 



(19) 



Under the previous approach we can determine the mean damage expansion rate rj 
which will be given by (here (...) are time averages): 



77 = {rj{t)) 

/ id(t+i : 

A |d(t) 
FXx(t))Qd(t) 
I d(t) I 



(20) 



(21) 



(22) 



such quantity can be easily computed by analyzing the statistical behavior of | 
F'(x(t)) d(t) |. This can be done by assuming that, on mean field grounds, 
F'(x(i)) can be replaced by a random matrix Q. The previous average (19) can 
then be estimated by considering the percent of I's in d{t) (i.e. | d{t) \ /N) and the 
same quantity for i + 1 (i.e. | d{t + 1) | /N). We have: 



d(t+l) 



1 - 



1 - 



2p{l-p)K I d{t) 



N 



N 



N 



(23) 



Now, in the thermodynamic limit {N oo) we get: 



r] — 1 — exp 



-2p{l-p)K^-^ 



N 



|d(t)| 
N 



2p{l-p)K 



(24) 
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This result could be derived in another way. By defining the normalized Hamming 
distance of a Boolean matrix fl as: 



N 



Ar2 51 



(25) 



where flij e {0, 1} We have: 



1 n d(t) 




n 1 


1 d(t) 1 




d{t) 








d{t) 





ft 1= 2p{l-p)K 



Prom (19) and (24), the Lyapunov exponent will be: 

A = \og\2p{l -p)K 



(26) 



(27) 



which determine the two classical regimes: A < (order) and A > (chaos) with 
the marginal case A = 0. In agreement with the boundary phase transition (2). 



6 Distance and Wolf's method 

This result has been consistent with the equation of distance evolution (1). If we 
interpret one of the replicas in the distance method as a perturbation of the another 
replica, the expansion in the time t of the perturbation will be: 

_ |rfi2(t + l) I _ 2p{l-p){l-[l-d{t)f} 

\d,,{t)\ - d{t) ^^^^ 

and approximating for small d{t): 

(1 - d{t))^ ^ 1 - Kd{t) (29) 

we have: 

ri{t) ^ 2p{l - p)K (30) 

i.e an expansion rate identically to (24). 

This approximation is prove sufficiently good as we show in fig. 5. The values of 
Lyapunov exponents through equation (27) with K = 3 and p variating and there 
equivalents values calculated through the distance equation (1) coincides. 

There are different methods to compute Lyapunov exponents in continuous sys- 
tems [4] . In order to show consistence we will demonstrate that it is possible compu- 
tate Lyapunov exponents from self-distance in consonance with the previous result. 
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The Wolf's method is used to numerically estimate Lyapunov exponents from 
time scries. In short, the method is as follow: get two points of the time series, let us 
say X(ti) and X(t2) and compute their relative distance: | X(t2) — X(ti) |. Assume 
that I X(t2) — X(ti) |< e, being e > very small. Next, compute the distance after 
T steps, i.e: | X(t2 + T) — X(ti +T) |. This time T is a fraction of the characteristic 
period or is defined in terms of the autocorrelation function. Repeating for n pairs 
of points and averaging, we obtain an estimation of the Lyapunov exponent: 

, 1 ,, jx(f; + r)-x(ti + r)| 

^ = S log |x(«-x(w (^1) 

For RBN, we can write an equation for the normalized Hamming distance between 
successive time steps in our system, i.e. the self-distance: dt^t-i- It easy to see that 
the self-distance is a new order parameter. This is a consequence of the combination 
of the distance method and the stable core. 
The iterated equation for the self-distance is: 

dt+i,t = 2p(l - - (1 - dt,t-i)''] (32) 

Wich formally is equivalent to equation (1) but different conceptually. The self- 
distance, as the stable core, does not require annealed replicas (as the distance 
method), and is computationally more easily to determine. In fig. 3-4 the numer- 
ical values of self-distance (points) are calculated in similar way that the distance 
(squares) with a very good agreement. 

If we approximate linearly close to the fixed point d* — 0, the function becomes: 

dt+i,t = 2p{l-p)Kdt,t-i (33) 
The iterated equation now is resoluble: 

dt+T,t = [M^-p)KYdt,t-i (34) 

Thus, we have: 



dt2+T — dt^+T 

dt2 ~ dt^ 



[K2p{l-p)f (35) 



i. e. aconstant value that, after introduced in the sum of (31), gives the Lyapunov 
exponents (27). 
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7 Summary 



In this paper wc have analyzed a new order parameter for RBN in terms of a f2- 
random matrix approach. Our order parameter deals with the percent of non-zero 
elements that is obtained from the hmit hm^^oo f2*. It is shown that the order 
parameter describes a (second order) phase transition at a critical point consistent 
with other previous analyses. 

An inmediate extension of the Boolean derivative approach is the construction of 
a measure for the damage expansion rate, rj(t) that give a quantitative characteriza- 
tion of how small perturbations propagate through the network. It has been shown 
that riit) provides a consistent measure of such sensitivity to spin flips. The time 
average over the Boolean products of the Jacobi matrix on the distance vectors can 
be successfully translated to a simple annealed method where only the statistical 
properties of the random matrix Q, matter. 

We also have calculated the Lyapunov exponents through the distance method. 
And we propose a new order parameter: the self-distance. This quantity opens the 
possibility of defining the Lyapunov exponent in a discrete system in analogy with 
the Wolf's method for continuous systems and the result is in agreement with the 
previously obtained. 
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